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Abstract 

In this paper we present a model independent analysis method following Bayesian statistics to search for 
neutrinos from point sources. The Bayesian approach allows us to obtain the full probability distribution 
function for both the background and the signal rate. As such, we have direct access to the signal significance 
and upper limits. Furthemore, it allows also to account for previous results via the concept of prior information 
without the need of the ad hoc application of trial factors. 

In this paper we derive an exact analytical expression for the probability density function for the source 
rate, and consequently also for the flux, of any steady point source assuming a uniform prior distribution. To 
investigate the validity of our Bayesian approach, we have applied this method to the public IceCube 40-string 
configuration data for 10 nearby blazars and have obtained a flux upper limit, which compares very well with 
the previously published result of IceCube, using the same data set. 

Key words: Neutrino Astronomy, Neutrino Telescopes, Active Galactic Nuclei, Gamma Ray Bursts, 
Bayesian Statistics. 



1 Introduction 

Bayesian approaches [1] are gaining more and more popularity in scientific analyses [2,3,4,5]. In this 
paper we outline a formalism following Bayesian inference to analyse signals from point sources, such as 
Active Galactic Nuclei (AGN) and Gamma Ray Bursts (GRBs), using data from neutrino telescopes. 

A frcqucntist approach is based on the long run relative frequency of occurrences in identical repeats 
of an experiment. Consequently, this can only provide the probability for a certain outcome under the 
assumption of a specific hypothesis. On the other hand, the Bayesian approach allows us to directly 
compute the probability of any particular hypothesis or parameter value based on observations. One of 
the great strengths of Bayesian inference is the ability to incorporate relevant prior information in the 
analysis. This provides a mechanism for a statistical learning process that automatically takes previous 
results into account. 

2 Bayesian Formalism 

We denote by p(A, B\C) the probability for hypothesis A and B to be true under the condition that C 
is true. The product rule [1], p(A,B\C) = p{A\C)p{B\A, C) = p{B\C)p(A\B,C), directly yields Bayes' 
theorem [6] : 

p{A\C)p{B\A,C) 

p(j4|s ' c) = — mc] — • (1) 

Bayes' theorem is extremely powerful in hypothesis testing. Consider a hypothesis H, some observed 
data D and prior information /. Bayes' theorem can then be rewritten as : 

p(D\I) 
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where 



p(H\D, I) = Posterior probability of hypothesis H. 

p(H\I) = Prior probability of hypothesis H. 

p(D\H,I) = Likelihood function, £(H). 

p(D\I) = Normalization factor. 

From Eq. (2) it is seen that the Bayesian formalism automatically provides a learning process. The 
first step is to encode our prior state of knowledge into a prior probability p(H\I). Based on the arrival 
of new data this is then converted into a posterior probability p(H\D, I). 

We will apply this approach to the analysis of astrophysical point sources. This method will enable us 
to obtain the full Probability Density Function (PDF) for the source rate from which we can derive the 
corresponding flux or an upper limit. 



3 Bayesian Assessment of Significance 

As outlined above, Bayesian logic allows us to make statements about the probability of various 
hypotheses in the light of obtained data. Following the developments described in [7], we quantify the 
degree to which data support a certain hypothesis and as such make an assessment of the significance. 
To quantify our degree of belief in a certain hypothesis H , we use the so called evidence 



e(H\DJ) = 10 log 



10 



P(H\D,I) 
p(H\D,I) 



(3) 



where H indicates hypothesis H to be not true. Due to the log 10 , this evidence reflects a decibel scale. 
In this context, the Bayesian observable t/j = — 10 log 10 p(D\H, I) and ip = — 10 \og 10 p(D\H, I) arc 
introduced. By combination of Eqs. (2) and (3) we obtain an expression for the evidence e(H\D,I) for 
any alternative to hypothesis H based on the data D and prior information I. In terms of the variables 
ip and ip, this yields : 

e(H\D, I) = e(H\T) + V - ? < e(H\I) + ip. (4) 

Eq. (4) implies that there is no alternative to a certain hypothesis H which can be supported by the 
data D by more than i[> decibel relative to H . 

Let us now consider an experiment where the probabilities pk corresponding to the various outcomes 
Ak on successive trials are independent and stationary. Such experiments belong to the so called Bernoulli 
class B m [1]. The probability p(ni....n m \B m , I) of observing nk occurrences of each outcome Ak after n 
trials is therefore given by the multinomial distribution : 

nJ 

p(D\B m ,I) =p( ni ....n m \B m ,I) = — — -.p n 1 \...plT ■ (5) 

m\....n m \ 

In terms of the Bayesian observable ip we then obtain for each H £ B m : 

m 

i'B,„ = -10[log 10 n! + ^(n fc log 10 p fc - log 10 n fc !)], (6) 

k=l 

which is an exact expression! 

In case the data is represented in histogram form, the above implies that n is the total number of 
entries, m represents the number of bins, is the number of entries in bin k and pk is the probability 
for an entry to fall in bin k. Once the various pu are known, the -0-value corresponding to a certain 
observed distribution can easily be obtained using Eq. (6). 

In our investigation of AGN, a histogram will be populated according to the angular difference a 
of each neutrino arrival direction w.r.t. the source location in the sky. The data from a homogeneous 
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background will follow a sin(a) distribution, reflecting the solid angle effect. This obviously determines 
all pk values once the angular bin size is chosen. 

By performing off-source measurements and using Eq. (6), we obtain the distribution of the correspond- 
ing background ip-v&hies by randomizing these background data multiple times according to the solid 
angle distribution. By comparing the Rvalue of the on-source data with the background '(/'-distribution 
we obtain the P- value or significance of our measurement. We can claim a source detection if the signifi- 
cance of a deviation from the background hypothesis has a P- value smaller than 5.73 TO -7 , corresponding 
to a 5<T effect in case of a Gaussian distribution. 

In case the observation does not lead to a significant detection, we determine an upper limit on the 
signal strength. To obtain such an upper limit we have worked out an exact analytical expression following 
a Bayesian approach using a uniform prior. The details of this procedure are provided hereafter. 



4 Bayesian upper limit determination 



To obtain an upper limit for a possible source flux, we first have to determine the upper limit for the 
source rate. Using Bayesian inference we obtain the full posterior PDF for the source rate and from that 
we can derive the corresponding flux PDF, via the concept of effective explained hereafter. 

In any experiment where events are recorded at a constant rate, the PDF for the number of observed 
events is described by a Poisson distribution. Our approach to obtain an upper limit is similar to the one 
outlined in [8] , except that we have obtained an exact analytical expression without any approximation. 

The data consist of on-source and off-source measurements, where the off-source data consist only of 
background events and the on-source data are a mix of background and source (i.e. signal) events. We 
start with the off-source analysis and use subsequently the obtained information as prior information in 
the on-source analysis. This implies that the current study is completely data driven and as such is a 
model independent search for a possible source signal. 



4-1 Off-source measurements 



Consider the case that in an off-source measurement N q r background events have been recorded over 
a time interval T g with a constant background rate b. Using Eq.(2) we obtain the posterior background 
PDF by : 

, u „ n p(b\I)p(N oS \b,I) (7 . 
^ |7W)= P (N oS \I) • (7) 

In the above equation the likelihood function p(N s\b, I) is given by the Poisson distribution corre- 
sponding to the measurement of N a g background events over a time span T q q at a constant rate b: 

(bT r F ) N °s e - bT ° 11 

p(N oS \b,I)= [bloS) e - . (8) 

Since the integrated PDF amounts to I (i.e. J b bmax p(b\N g, I) db = I), the normalisation factor p(N Q g\ I) 
appearing in Eq. (7) is given by : 

p(N oS \I)= [ ^ p(b\I)p(N oS \b,I)db. (9) 

Jb min 

As mentioned before, p(b\I) is the prior PDF for the background rate. For our analysis we use a uniform 
prior [8], which is given by 

PW) = 7 ^j— • (10) 

"max "mm 

The uniform prior attributes the same probability to each value of the background rate within the 
indicated range, reflecting our total ignorance about the actual background rate. 

To cover the full range of possible background rates, the minimum value of the rate b is taken to be 
zero and Eq. (10) can be written as 

p(b\I) = J-. (II) 
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Using this expression for p(b\I) together with Eqs. (8) and (9) we obtain an analytical expression for 
the normalisation factor 

/■bmax i ( bT \N o[[ -bT o[[ 

p(N on \I) = / ^- ^r; f db. (12) 

JO "max JV ff! 

Solving the above equation (for details see Appendix Al) we get : 

p{N oS \I) = — — , (13) 

Omax i Voff ! J- off 

where 7(0, x) is the Incomplete Gamma function [9]. Substitution of Eqs. (8), (11) and (13) in Eq. (7) 
yields 

p(b\N oS ,I) = — -, (14) 

r Y( N otf + l,O m axToff) 

which represents the posterior background rate PDF. 

4-2 On-source measurements 

Consider the case that in an on-source measurement N on events, consisting of signal and background, 
have been recorded over a time interval T on with a constant signal rate s and background rate b. Following 
Eq. (2) the joint probability of source and background is given by : 

p{s, b\N on , I) = — . (15 

P{Non\I) 

Using the product rule [1] we can write the above equation as, 

/ ti at- n p(b\I)p(s\b,I)p(N on \s,b,I) 

p(B,b\Non,I) = , (16) 

where p{b\I) is the prior probability for the background rate, which is in our case the posterior back- 
ground PDF obtained from the off-source measurement reflected in Eq. (14). The likelihood function 
p(N on \s, b, I) is the Poisson distribution for the combined signal and background rate (s + b). The nor- 
malisation constant p(iV on |/) is obtained, as outlined in the previous section, by integrating the numerator 
of Eq. (16). 

It is important to note that since the source rate s and the background rate b are independent, we can 
write p(s\b, I) = p(s\I). Like for the background case discussed before, we use a uniform prior for p[s\I), 
i.e. 

p(s\I) = — =p(s\b,I). (17) 

° max 

As before, this reflects our total ignorance about the source rate. 

By substituting in Eq. (16) the expressions of Eq. (17) for p(s\b,I), (14) for p(b\I) and (8) for a total 
rate (s+b) we obtain the joint PDF for the source and background rates : 

jJ^CCU"* e ~ bT °" ■< i + »)""" e-' b * ,,T " did. 

However, we are interested in the posterior PDF for the source rate alone, independent of the back- 
ground. The Bayesian formalism allows us to obtain this posterior PDF by marginalisation [1] of the 
joint PDF, Eq. (18), with respect to the background i.e. 

p{s\N on , I) = [ "" X p(s, b\N on , T)db. (19) 
Jo 

Solving the above integral we obtain an exact expression for the source rate posterior PDF (for details 
see Appendix A2) : 

sTon V^^on S i (r op +r off ) i 7 (iV-»+l,Mm, 1 x) 

/ I at n _ C f-j=0 i!(AT on -i)! 

P{S\l\ on , 1) - Nm (T on +T off )j7(Af-j + l. Mmax ) 7 (j + l, Smax T on ) ' ^1 
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where u max = 6 max (T on + T oS ) and N = N on + N otf . 

In case no significant source signal is observed, Eq. (20) allows us to derive an upper limit for the 
source rate. As an example, the 90% source rate upper limit s u .i. is given by : 

Ul 'p(s|iV on ,/)ds = 0.9 (21) 



5 Analysis of ten nearby blazars 

The approach presented in the current paper is tailored for sources with a constant signal rate. As 
such, our primary goal is the analysis of AGN which are expected to be steady sources of high-energy 
neutrinos and of which the locations are well known. By selecting a small region around each source 
location, several sky patches are defined from which data were collected to search for a possible deviation 
from the background "noise" . The background is determined by performing observations of similar sky 
patches which do not coincide with AGN locations. 

To validate the analysis procedure described in this report, we use the public data [10] of the muon 
neutrino candidate events recorded by the IceCube Neutrino Observatory [11,12] in its 40 string con- 
figuration (IC40), that collected data during the season 2008-2009. Our analysis was performed on ten 
nearby blazars (a special class of AGN with one of the jets pointing in the direction of the Earth) follow- 
ing the approach described in the previous sections. The blazars were selected from the online "Roma 
BZCAT Multi-frequency Catalogue of Blazars" [13] and are listed in Table 1. These blazars were chosen 
to be nearby, i.e. with a small redshift, and in the Northern hemisphere, to reduce the atmospheric muon 
background for the IceCube measurements. 

We investigated each source individually and selected the recorded events within a 5° window around 
the source position. The 5° search area implies selecting all events within the 5<r angular uncertainty of 
the Icecube track reconstruction, which was determined to be of the order of 1° [12,14]. Note that the 
uncertainty on the blazar positions is less than an arcsecond. 



Blazar Name 


Right Ascension 


Declination 


Redshift 


Distance [Mpc] 


BZUJ1148+5924 


177.20983 


59.41567 


0.011 


46.2 


BZUJ0048+3157 


12.19642 


31.95697 


0.015 


63 


BZUJ0319+4130 


49.95067 


41.51169 


0.018 


75.6 


BZUJ0319+4130 


107.39246 


50.18225 


0.02 


84 


BZUJ0153+7115 


28.35771 


71.25181 


0.022 


92.4 


BZUJ1719+4858 


259.81025 


48.98042 


0.024 


100.8 


BZUJ1632+8232 


248.13321 


82.53789 


0.025 


105 


BZUJ1715+5724 


258.84579 


57.41119 


0.027 


113.4 


BZBJ1104+3812 


268.95183 


62.61225 


0.027 


113.4 


BZUJ1755+6236 


166.11379 


38.20883 


0.03 


126 



Table 1. Nearby blazars [13] used in the current analysis. 



5.1 Bayesian assessment of significance 

As outlined in [7], stacking the recorded events within the 5° cones of these ten blazars, according to 
their angular distance a from the actual blazar position, we obtain the distribution of events presented 
in Fig. 1. As mentioned in the previous sections, we need an expression for the various probabilities pk of 
Eq. (6) to derive the -0-value of the signal. In our case, the background is homogeneous and consequently 
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the probabilities have to be consistent with the solid angle effect within our 5° cone, i.e. 



Pk = 



1 - cos(5° 



[cos(w • k) — cos(w ■ (k + 1))] , 



(22) 



where w is the width of each bin in our stacked histogram. 

The stacked data, reflected in Fig. 1, comprised 407 events recorded over a time period of 375.5 days 
[15]. We have chosen 0.1° as angular bin width, which is well below the IceCube angular resolution of 
1° so that our analysis can be regarded as unbinned. Using Eqs. (6) and (22) for a number of entries of 
n = 407 and a number of bins m = 50, the data represented in Fig. 1 yield ^observed = 456.57 dB. 



o 
u 
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Fig. 1. The stacked distribution of events within a 5° cone of all our 10 Blazars of Table 1, with a the angle 
between the corresponding blazar location and the reconstructed arrival direction. 
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Fig. 2. The stacked distribution of events within a 5° off-source background region, 180° shifted in azimuth from 
the various blazar positions. 

To quantify our belief that the stacked histogram represents only background, we need to compare this 
Vvbserued-value to the distribution obtained from the background only data. As the efficiency of IceCube 
is zenith angle dependent, we performed measurements of the background rate in 5° regions of the sky 
shifted from the various blazars positions only in the azimuth angle, keeping the zenith angle constant. 
However, the specific IC-40 configuration of IceCube is also azimuth angle dependent, so we made shifts 
of 180° in the azimuth angle to eliminate the azimuth angle dependence. The IC-40 sample has been 
taken over a detector live time period of 375.5 days, so that both the exposures for on-source, T on , and 







off-source, T Q s, amount to 375.5 days. The off-source measurements yielded a total of 409 events. These 
stacked background measurements are shown in Fig. 2. 

Randomising the background data 10 6 times according to the solid angle distribution, we obtain a 
distribution for the background tpbkg -values, which is presented in Fig. 3. 
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Fig. 3. The distribution of ^-values of 10 6 randomizations of the background measurements. 

Comparison of the actual observation ^observed with the background distribution ipbkg yields the P- 
value of our measurement, which amounts to a P-value of 0.87. This implies that we cannot claim the 
presence of a signal. Consequently we will proceed to give an upper limit on the signal strength. 

5.2 Upper limit determination using Uniform priors for Source and Background 

5.2.1 Determination of the background rate 

The posterior background rate PDF is obtained by substitution of the previously mentioned values 
of T ff and N Q g_ in Eq. (14) and by using a sufficiently large value 6 max = jr^- • 100 = 1.16 mHz. The 
resulting background rate PDF is shown in Fig. 4. 
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Fig. 4. Background rate PDF using a Uniform prior. 

5.2.2 Determination of the source rate 

The posterior source rate PDF is obtained by inserting the previously mentioned values of T on , T Q ff, 
A^on and N g in Eq. (20) and by using a sufficiently large value s max = 1 Hz. The resulting source rate 
PDF is shown in Fig. 5. 
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Fig. 5. Source rate PDF using a Uniform prior. 

Using the PDF shown in Fig. 5 and applying Eq. (21), we obtain the 90% upper limit for the source 
rate : 

s u .i. = 1.41 x l(T 6 Hz. (23) 

Note that this rate does not take into account the reconstruction efficiency. The latter is taken into 
account by converting the source rate upper limit into a flux upper limit by means of the so called 
Effective Area, A c ff , which is defined as 

number of observed events 

Acff = ; ■ (24) 

incoming flux 

By definition A e ff takes all reconstruction efficiencies into account. For the current analysis we use 
the angle averaged Effective Area, weighted by the observed energy estimate [10], which corresponds to 
A ff = 9.95 x 10 5 cm 2 . Our analysis was performed on a circular area of 5° centered on each of the 10 
sources, representing in total 10 x 0.0239 sr= 0.239 sr. Together with the expressions (23) and (24) this 

leads to a 90% upper limit for the signal flux of $ u .i. = — = 5.93 x 10~ 12 s" 1 cm~ 2 sr _1 . 

A e ff • 0.239 

However, this flux upper limit does not take into account the effect of neutrino oscillations. At the 
source we have : v e : v T = 2 : 1 : [15]. Assuming maximum oscillation we will observe on the Earth 
i/p : v e : v r = 1 : 1 : 1. A tiny fraction of the v T will produce a muon which might also be detected 
in IceCube and as such have entered our event sample. However, we will neglect this effect since it is 
marginal and would require a special simulation which is beyond the scope of this paper. So our final 
value for the 90% upper limit for a signal flux is : 

$ u .i. = 11.86 x 10~ 12 s' 1 cm- 2 st- 1 . (25) 

This upper limit may be directly compared to the result published by the IceCube Collaboration [15], 
which was obtained with a different analysis concerning a search for a diffuse high energy neutrino flux 
in the full Northern hemisphere: $ u .i. = 8.89 x 10~ 12 s _1 cm~ 2 sr _1 . 

The limit obtained in this paper compares rather well to the above, and therefore indicates that the 
analysis method presented here is suitable to investigate neutrino fluxes from point sources. 



6 Conclusion and Outlook 

In this paper we have presented a mechanism to analyse point sources using data from neutrino 
telescopes following Bayesian inference. We have shown how to test a background only hypothesis using 
Bayesian statistics for sources like GRBs and AGN and we have provided a method to determine a 
possible constant source rate. For the latter we have worked out an exact analytical expression using 
uniform priors, that can readily be used to obtain the probability density function for the rate of any 
steady source. We have also shown how to obtain upper limits for the corresponding flux in case the 
observation does not lead to a significant signal. 
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Applying this method we analysed the public IceCube 40-string configuration data for 10 nearby 
blazars located in the Northern sky. We have found that the on-source data is consistent with a back- 
ground only hypothesis. Therefore we have determined a 90% upper limit. Our upper limit agrees well 
with the previously published IceCube upper limit based on the same data set. 

Furthermore, it is our intention to extend the current method also for non-steady sources like for 
instance GRBs. Apart from providing a signal significance for discovery [7], this should also provide a 
mechanism to accurately determine flux upper limits. The latter requires a non trivial extension of the 
method presented here, which is currently under study. 
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Appendix Al 



The normalisation factor of the posterior background rate PDF given in Eq. (12) may be written as : 

P( N ° S l J ) = I / ( bT os) N ° S e ~ bT ° S db - ( 26 ) 

Omax ^ S- Jo 

The integral part can be expressed as the so-called Incomplete Gamma function given by : 



7(0,1) = / e-H°- l At. (27) 
Jo 

Using this expression, we can rewrite Eq. (26) as follows : 

^ off |/) = ^Ar3 / bmaxT °^^ off )^ e --o ff ^k) (28) 

(29) 



frmax N oS ! J T oS 
1 l(N oS +l,b ma , x T oS ) 



&max N o^- T oS 
which is the expression reflected in Eq. (13). 

Appendix A2 

According to Eq. (19) the posterior PDF for the source rate alone is given by : 

rbmax 

p(s\N OD ,I) = p(s,b\N on ,I)db, (30) 



Jo 

where 

/ MKr n P (s\b,i)p(b\i)p(N on \s,b,i) 

p{s,b\N on ,I) = — . 

P{Non\I) 

Substitution of the various expressions given in section 4 yields : 

1 T off (6T off ) JV offe-'' r off ( b+s )«on To " n on e -( &+s )T on 



p(8, b\N on , I) = SM 7(^V off +l, bm axT off ) W _ (31) 

rb max rs max 1 T oS (bT oS ) N otte-" T oti (H^^C^-tH-lTon 
JO JO smax ' j(N oS +l,b max T oS ) ' N on i Q0 QS 

Combination of Eqs. (30) and (31) yields : 

, 6max _JL_ . T nff (bT nff foSe- bT oS (6+a) iVon r< ffin e -( t+ .)T on 

v(s\N I)— smax 7(JV off +l,h m axT off ) jVpn! ^ 

rfemax f s m ax 1 T oS (bT oS ) N oSe-" T oS (b+s)~onT f n '"e-(''+°>?bn 
JO Jo smax ' 7(A r ff+l,f'maxT off ) ' JVon! 







rfcmax _j_ . T off (bT off ) Jv off e ^"off _ (H^^C^-^'^ dfo 

Jo smax 7(JV of j+l,hmaxT off ) Won! 



rbmax fs m ax 1 T off (br off ) W off e'" T off (b+s)~onT^ n " e -("+°>?bn 

Jo Jo s max ' y(N oS +l,b max T oS ) ' JV on ! 

A 
B 



(32) 
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Considering the numerator of Eq. (32), we obtain 



_ J_ T off (6r off ) jV offe- &T ofF (b + 8 )^onr o iVon e -(6+ a )ron 

JO "max 7 (A/off + 1 , 6 max T of f ) iV on ! 



1 rpN nT N ° s+1 1 f'max 

1 J on J off ' 



N on \ 7(-^off + 1, kmaxT'off) J 



~(b + S ) N ™b N offe-( b+ ^ T ™e- bT offdb. 



Using the Newtonian Binomial, (a + b) n — X)"=o i\(n-i)\ a " we nn< ^ ■ 

i 7^A'onT Aroff+1 1 /V I o» p-sTon /•'max 
A = — — io " V °" / h Non+N oS -i -b{Ton+T oS ) d& 



A^on! 7(^off + 1, b max T oS ) f^ Q i\(N on - i)\ J 
Using the Incomplete Gamma function, Eq. (27), and simplifying the above equation, we finally find 

1 T JVon7i Ar off+ 1 ^on j - s T on ~( AT _ A i 1 „ \ 

^ _ f on J off s e 7U V z + MWxj /gq\ 

~ Smax 7(^off + 1, &maxT off ) ^ i\(N on - »)! (T on + T a e) N ^+ N oS-*+^ ' 1 j 

where w max = o max (T on + T off ) and iV = iV on + N oS . 

Applying the same procedure, we obtain for the denominator : 

/•brnax /-smax 1 T oS (bT oS ) N ott e - &T oft (6 + s )^on I W one -(&+ s )T on 

B = / / ■ — — — • ; do d.s 

Jo Jo s max y(N ff + l,b 

ma. x Toff) AU! 

s max 

ATds 



i 7iA'on7^ 7V off+ 1 ™2? 1 -vflV-i + l i/ "| f s max 

Smax 7(A'off + 1, OmaxT off ) ^ j!(iV on - j)l (T on + T oS ) N ™+ N o«~^ J Q 

1 T N n °»T^ +1 g 1 7(JV-j + l,w) 7(j + l,s max T on ) 



.s 



Smax 7(A r off + 1, OmaxToff) ^ j!(iV on - j)l (T on + T off )ATon+iV off -j + l ^ + 1 

Which finally gives, after simplification, for the posterior source rate probability density function 

-sT on X^Non s' {T on +T oS Y l(N-i+l,u mSLX ) 
\pivon (Ton^-i ff; J Tuti,smaj(Jo]i; T 



p - = j- on V on — 

M/v n= B i\(N on -i)\ , . 

P{S\l\ on ,l ) „Af on (T on +r off )j 7(3 + 1, smax Ton) 7(^-3 + 1, «max) " ^ > 
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